Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

make unitcell matrix column major #71

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

make unitcell matrix column major #71

wants to merge 1 commit into from

Conversation

lmiq
Copy link
Contributor

@lmiq lmiq commented Jan 26, 2023 •
edited
Loading

The docs (here) imply that the matrix returned by matrix is column major. However, this is not what is actually happening. (example file attached).

Example:

julia> traj = Chemfiles.Trajectory("./traj_comp.xtc");

julia> frame = Chemfiles.read_step(traj,0);

julia> unitcell = Chemfiles.UnitCell(frame);

julia> Chemfiles.matrix(unitcell)
3×3 Matrix{Float64}:
 77.395     0.0     0.0
 38.6975   67.026   0.0
 38.6975  -22.342  63.1928

as such, this matrix does not satisfy the conditions gromacs requires.

Just transposing the matrix would be easy to solve as Matrix(transpose(Chemfiles.matrix(...))). However, I'm not sure if it makes the definition inconsistent with the return of the unit cell in terms of lengths and angles. For instance, in this trajectory, we get:

julia> Chemfiles.matrix(unitcell)
3×3 Matrix{Float64}:
 77.395     0.0     0.0
 38.6975   67.026   0.0
 38.6975  -22.342  63.1928

julia> Chemfiles.lengths(unitcell)
3-element Vector{Float64}:
 77.39500045776367
 77.39500188456586
 77.39500082545612

julia> Chemfiles.angles(unitcell)
3-element Vector{Float64}:
 90.00000019466002
 60.000000157157025
 60.000000609835716

But if we tranpose the matrix, we do not recover these angles and lengths:

julia> mt = Matrix(transpose(Chemfiles.matrix(unitcell)));

julia> ut = Chemfiles.UnitCell(mt)
Chemfiles.UnitCell(Chemfiles.CxxPointer{Chemfiles.lib.CHFL_CELL}(Ptr{Chemfiles.lib.CHFL_CELL} @0x000000000910c100, false))

julia> Chemfiles.angles(ut)
3-element Vector{Float64}:
 108.4349482115009
  65.90515744788931
  75.03678223987454

julia> Chemfiles.lengths(ut)
3-element Vector{Float64}:
 94.78912988199575
 70.65164746673275
 63.19275379180908

The simplest alternative would be to fix the docs and indicate that the matrix is row-major. That is somewhat inconsistent with Julia, but then at least it lets to the user the decision to transpose the matrix and arch with the consequences.

The other alternative is to check if this PR is fine, I have only transposed the matrix back and forth when converting it from the UnitCell.

traj_comp.xtc.zip

@Luthaf
Copy link
Member

Luthaf commented Jan 29, 2023

Right, this is an issue that comes from the C++ side of things. Over there, the data is stored as row major, but containing the cell vectors as columns (which is confusing, and which I want to change in the next breaking release).

When passing data to Julia, the data is interpreted as column major, meaning the cell vectors are now rows of the matrix ...

I think this PR is a nice way of handling this disparity until we update the representation on the C++ side. I'm not sure if I should put these changes into a new release though, since users of Chemfiles.jl might already be handling the rows/column disparity, and this would be a breaking change for them. So you would need to pull in chemfiles from git to have these changes.

The alternative would be updating the docs for now, and waiting for the next breaking release to change the behavior. Do you have a preference about this?

@lmiq
Copy link
Contributor Author

lmiq commented Jan 29, 2023

I thought about that. It is not clear to me if that can be considered a bug and the docs are correct, though.

I don't have a strong opinion. Feel free to deal with it in the way you find safer.

@lmiq
Copy link
Contributor Author

lmiq commented Feb 14, 2023

I sort of need a definition on this. Lets keep it as is, and let the change in row order for the next breaking version? Then I'll handle that within my package.

@Luthaf
Copy link
Member

Luthaf commented Feb 14, 2023

That works for me. Sorry, I should have told you I'm on holidays right now and wanted to address this after I'm back.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers
No reviews
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.
2 participants
@lmiq @Luthaf